Electrical field induced shift of the Mott Metal-Insulator transition in thin films 
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The ground state properties of a paramagnetic Mott insulator are investigated in the presence 
of an external electrical field using the inhomogeneous Gutzwiller approximation for a single band 
Hubbard model in a slab geometry. The metal insulator transition is shifted towards higher Hubbard 
repulsions by applying an electric field perpendicular to the slab. The spatial distribution of site 
dependent quasiparticle weight shows that the quasiparticle weight is maximum in few layers beneath 
the surface. Moreover only at higher Hubbard repulsion, larger than the bulk critical U, the electric 
field will be totally screened only for centeral cites. Our results show that by presence of an electric 
field perpendicular to a thin film made of a strongly correlated material, states near the surface will 
remain metallic while the bulk becomes insulating after some critical U. In contrast, in the absence 
of the electric field the surface becomes insulating before the bulk. 

PACS numbers: 71.30.+h, 71.27.+a, 73.61.-r 



I. INTRODUCTION 

The rich physics of strongly correlated materials in 
combination with the need to overcome the scaling lim- 
its of current silicon based semi-conductor materials in 
microelectronic industry has resulted in an increased ac- 
tivity in this field. Special attention has been focused on 
vanadium dioxide(y02) which shows an abrupt Metal 
Insulator Transition(MIT) near room temperature due 
to a structural phase transition.— One has found that an 
electrical field is able to trigger MIT in VO2 , without any 
structural transition, which is mostly dominated by elec- 
tron correlations rather than a Peierles distortion^. Also, 
a first order MIT is observed by applying an electrical 
field in a two terminal model of VO'^^ without electri- 
cal breakdown of the material. Note that an, electric 
field driven MIT and metal-superconductor transitions 
have been observed at the interface between LaAlOj, and 
SrTiOs.— These kind of transitions may be related to 
a charge transfer mechanism. A nonlinear dependence 
of the conductivity on the electrical field is reported 
for the highly correlated transition metal chalcogenide 
Ni{S, 5*6)2 ^^'^ ^ continuous MIT is observed in this 
case^. 

In this paper we investigate the behavior of the ground 
state of a single band Hubbard model^ in the presence 
of a perpendicular electric field by using the Gutzwiller 
approximation (GA)^. Originally GA is rooted in the 
Gutzwiller wave function used to reduce the contribution 
of high energy states due to Hubbard repulsion and it was 
shown to be exact in the limit of infinite dimensions!^"— 

While an analytical solution exists only for one 
dimensiortii, in comparison to other approximate meth- 
ods the GA is equivalent to a slave boson mean field 
approach(SBMF)^ for zero temperature but in contrast 
to dynamical mean field theory (DMFT>i^, it is not able 
to give any information about higher and lower Hub- 
bard bands. Instead it gives a reasonable understanding 
about the low energy excitations near the Fermi surface^ 
by supplying the quasiparticle weight of electrons such 
that one is then able to describe the mobility of elec- 



trons. Also, GA cannot give any information about the 
insulating state, instead we are only able to investigate 
the properties of the system by approaching the tran- 
sition point, Uc, from below^^. This method was used 
by Brinkman-Ricei^ to investigate the MIT of the sin- 
gle band Hubbard model and it allowed them to predict 
the critical Hubbard repulsion which is finite in two and 
three dimensions ([/c = 16t for 3D). While not as accurate 
as DMFT, GA is less computationally intensive and thus 
allows the description of inhomogeneous systems such as 
thin films subjected to a perpendicular electric field. 

Although our simplified approach is only qualitative, 
it gives important information about how one may be 
able to spatially tune the quasiparticle weight distribu- 
tion near surfaces and interfaces. This could be rele- 
vant for future studies; for example for an inhomoge- 
neous bad metal-superconductor transition by the charge 
transfer mechanism which may be responsible for the SC- 
Insulator transition observed at the interface of a band in- 
sulator and a strongly correlated material^. We will show 
that by applying a perpendicular electric field, charges 
will be trapped at the surface of the Mott insulator and 
shift the MIT for the surface states. 

The outline of the paper is as follows. In section II 
we review the concept of GA and how the inclusion of 
on site potentials may change the situation. In section 
III we introduce our model for the slab geometry, present 
the numerical scheme used and analyze the corresponding 
results. Finally in section IV we present our conclusions. 



II. GUTZWILLER APPROXIMATION IN THE 
PRESENCE OF AN ELECTRIC FIELD 



In order to address the narrow band efi'ects in transi- 
tion metals with d or f orbitals for which correlation ef- 
fects play a major role in the behavior of the system the 
simplest model that is able to explain the most impor- 
tant terms of the Coulomb interaction between electrons 
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is the well known Hubbard model, 



+ C^j^Cia) +'YUflianis. (1) 



<ij><T 



We will describe the ground state properties of the Hub- 
bard model by using the Gutzwiller approximation which 
suppresses the contribution of high energy configurations 
(here configurations with higher number of double occu- 
pancies). This is done by introducing a trail wave func- 
tion which contains variational parameters to be used 
subsequently to minimize the total energy of the system. 
Our aim is to investigate the properties of a strongly 
correlated system in the presence of an external electri- 
cal field which will appear in the Hamiltonian as a posi- 
tion dependent potential. The induction of such an inho- 
mogeneity is not random and we still have translational 
invariance in the direction perpendicular to the applied 
field. To study the ground state properties in the ab- 
sence of the electric field, the Gutzwiller wave function is 
defined as: 



IV^o), (2) 



where the double occupancy operator is Di = fii„hi, 



the 



variational parameters Qi are introduced to reduce the 
contribution of high energy configuration's in the many 
body wave function \^g), and is the unprojected 
non-interacting (Fermi sea) many body wave function. 
Although it is obvious that by the inclusion of on site 
potentials no new variational parameters are needed be- 
cause they do not induce any new correlations since the 
term is a single body interaction, nevertheless we will 
prove it rigorously. To obtain the normalization factors 
in the limit of spatial infinite dimensions, for which the 
Gutzwiller approximation is exact^i^^, we have to remove 
spatial correlations which occur in infinite dimensions to- 
gether with on-site Hartree contributions which remain 
in the d = oo limit. This can be done by introducing 
a new expansion parameter following the guidelines of 
Ref [l^. If we include on-site potentials for capturing 
the effects of external fields, the Hamiltonian becomes: 



E 



(3) 



In order to find the ground state of the Hamiltonian in 
Eq. ([3|) we introduce new variational parameters, Q„ and 
Cig. , to decrease the weight of the occupancy of the sites 
with higher on-site energy. The Gutzwiller wave function 
now becomes: 



llpg) = [1 - (1 - Qa)n.,a] [1 - (1 - C^a)n^a] X 



f-(l-.g,)A Ha) 



(4) 



The standard way of removing on-site Hartree contribu- 
tions is to introduce the fugacity factors and ^J-ig^, 
the expansion parameter Xi and the non interacting state 



\(po)- Then the Gutzwiller wave-function can be written 
as: 



'Pa) 



Y[{l+x,iD,-D"^))\ipo). 



(5) 



The Hartree double occupancy operator can be defined 

- HF 

as A = h^a{nis)Q + {nicr)onig - (ri,cr)o(^t<j)o and it is 
the result of the usual mean field decomposition hi^r — >■ 
fiia- - (^^^o■)o• By defining Qa- = g/'" , (is- = 9t^"° ^ M^cr' = 
Act + fJ-ia and fii^' = /3,g + fi^g- we have: 



(6) 



Therefore by using the above change of variables it is pos- 
sible to obtain the same renormalization factors for the 
infinite dimensions limit as stated iniS. Moreover by us- 



ing the condition {fiicr) 



which holds for infinite 



dimensions it can be inferred that the physical counter- 
parts of the new variational parameters, Q^, are {nia-)Q, 
In minimization procedure we need to minimize the en- 
ergy with respect to \ipo) together with local variational 
parameters gi that one needed to describe the correlation 
effects. In short, the addition of on-site potentials does 
not add any new variational parameters and the proce- 
dure of finding the ground state is the same as in the 
conventional Gutzwiller method. Thus the expectation 
value of the Hamiltonian: 



(7) 



has to be minimized only with respect to gi and \ipo)- 
Here the renormalization factors qi„ depend on the local 
density of the non-interacting state \ipq) and gf. 



1 



Qicr 



("io-)o(l - {nia)a) 



V di{{hia)o - di) + \J {{niij)o - di){l - n,^o + dj) 



(8) 



where n^^o = ("■1(7)0 + (^iff)o' '^hil^ gi are described by the 
following equations which holds in infinite dimensions: 



9i 



dijl - Tiifi + di) 
i{nia)o - di){{nia)o - di) 



(9) 



Although in normal metals any deviation from half fill- 
ing may lead to the lowering of the electron conductiv- 
ity, in strongly correlated materials these deviation play 
different role because of the dependence of renormaliza- 
tion factors of tight binding parameter, 9^0-, on the local 
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charge density. Thus one may predict that if an applied 
electrical field would be able to change the charge distri- 
bution of the system then it will be able to change the 
electron conductivity and even shift the metal-insulator 
transition point. 

In practice minimizing the expectation value of the 
Hamiltonian is difficult because of the existence of a large 
number of variational parameters in \ifQ) together with 
the dependence of the renormalization factors on 
This will lead to a highly nonlinear set of equations. In 
order to alleviate some of the difficulties it is possible to 
allow local densities and |(/5o) to vary independently in 
the minimization procedure. Then by introducing \i„ as 
Lagrange multipliers it is possible to ensured that the lo- 
cal charge densities of the Gutzwiller wave function are 
equal to the local charge densities of the non interact- 
ing state. Other multipliers. A and are introduced in 
order to ensure total charge conservation and guarantee 
that |(po) is normalized. Therefore the final form of the 
energy expectation value is: 



~A{N -Y,na^) + E{1 - i^ol^o)), 



(10) 



where t, 



\/Qia\/Qja-tij are the renormalized hopping 
amplitudes. To find the optimum energy first of all we 
vary the \ipo) for which we have the following Schrodinger 
like equation: 



-h.c.)\Lpn)+22,{'^i+^i<y)f^i'y\Vo) = E\tpo), 

i 

(11) 

which has to be diagonalized for both spins. The sums 
are up to the filling of the system. This non-interacting 
energy is the amount of kinetic energy which is stored in 
the quasiparticle state |<po)- Then \ipo) is substituted in 
ea.lfn?)) and the expectation value becomes: 



{H) = ENi + J2ud, + AiN-J2n.^a) + J2^^ 



(12) 



where Ejsij is the non-interacting energy which depends 
on the variational parameters n.i^, \i„ and \^po). |<po) is 
now a function of the variational parameters \i„ , Uic, and 
gi, and the above energy functional has to be minimized 
in accordance to all these parameters. This leads to the 



following set of saddle point conditions: 

d{H) 



dk 

d{H) 

d{H) 
d{H) 



0, 
0, 
0, 
0. 



(13) 



III. MODEL AND NUMERICAL SCHEME 
A. Model 

Our model is a slab geometry in which we have trans- 
lational invariance in x and y direction and finite size in z 
direction. In addition we apply a linear potential profile 
from —v/2 to +v/2, in the z-direction. With the above 
assumptions the expectation value of the Hamiltonian 
can be written as: 



{H)^ {(fiol ^ {~2tq^a{cosk^ + cosky) + Vi + \ic,)c\^^^ 

Z,A;|| ,(T 

- 51 V^v/^*(cL||aCjfe||<T +4fc||<TC»fc||'T)|V50) 

<ij>k,,a 



N) + Eil- < vjol'/'o >) 
(14) 



where i and j correspond to atoms in the z direction and 
Nk^^ — Nk^Nky is the total number of k points. 

First we minimize the energy with respect to \(po) 
which leads to the following eigenvalue problem: 

{-2tqi„{coskx + cosky) + Vi + \ia)c\,,^^a^ik^\'y\'^Q > ^ 

fc|| ,(T 

?J^^(4-|,<TCjfc||^ + c]^ha)Wo >= £^fe|| IVO > 



<ij>fe||Cr 



(15) 



Eq. has to be solved for each fc|| point and in order 
to find the non-interacting ground state the eigenvalues 
will be summed up to the desired filling level: 



E 



NI 



E ^'^ii'"' 

E<Ef 



(16) 



where Ep is the Fermi energy of the quasiparticle states 
and n is the quantum number for the energy level of each 
k point. 

In the next step the above non- interacting state |</Jo >j 
which is now an implicitly function of all variational pa- 
rameters Aio- , Ilia J Qi and should be inserted into 
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Eq. dl): 

(H) =< ip[Xia,n^cr, 9t, A]\Ho\ip[Xia-,ni^, gi, A] > 

isCT 2,(7 



which holds when the wave-function is an eigenfunction 
of the non-interacting Hamiltonian and obtain the fol- 
lowing set of saddle point equations for the paramagnetic 
(17) case (< hia- >=< hig- >): 



We therefore minimize the total energy according to the 
variational parameters by considering: 

^ < m\Ho[x]\m >=< m\Q^Ho[x]\m >, (is) 
I 



-5, 




(19) 



d<H> 



2 < >fo\^^~2t{cosk,^ + cosfc.y)^^c]'j.||^Cifc||<^ 



(20) 



d<H> 

dXia 



('^Ol E 4,|aC»fc||<T|Vo) - A^fc||rii,T = 0, 



(21) 



d<H> 
dA 



(22) 



In order to numerically solve the above set of non-linear 
equations, we use MinPack. l^^i-*^^ which uses a trust- 
region-dogleg method, while for the k-space summation 
we choose a 8 x 8 Monkhorst-Pack^^ k-grid for which the 
energy is well converged for this kind of grid. From the 
above equations it is obvious that the Jacobian matrix 
required by the nonlinear solver has to be calculated by 
a finite difference method because no analytical evalua- 
tion of the Jacobian matrix is possible. Also note that 
the Jacobian matrix is dense and all of its elements are 
nonzero. 

We also tried to implement another approach by solv- 
ing Eqs. (|15p and (|19p - (|22|) iteratively by starting from 
an estimation of the variational parameters and a cal- 
culation of |(/3o > which are then supplied to the set of 
Eqs. (Unil-imi) to find a new set of variational parameters 
and then repeat the whole procedure. The iterative ap- 
proach did not converge for values of J7 > 4< which could 
be because of the high non-linearity of the equations for 
large U. Other authors also reported similar problems 



r 



with such an iterative scheme^^. 



Although the second approach is less costly, because 
the Jacobin matrix in the first method is updated at each 
variation of the parameters it is more likely that the first 
method converges better particular for large U when we 
have a large dependence of |iy9o > on the variational pa- 
rameters. In the next sections we report results for qt 
as the position dependent quasiparticle weight which is 
an indication of the mobility of the electrons in Fermi 
lic^uid theory. It is possible to show that the inverse of 
this factor is proportional to the mass renormalization 
which is divergent for qi = and which corresponds to 
an insulating phase^^. The quantity Vi = vi + Xia -I- A is 
considered as an effective potential which acts effectively 
only on \^po). The parameters U and v are scaled with 
the tight binding parameter t. 
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FIG. 1. (a) Quasiparticle weight distribution for U < 16t, 
7V2 = 100 and v = 2t and v = 0; (b) charge distribution for 
U < 16t, iV^ = 100 and v = 2t. Note that for d = the 
system is at half-filling, = 1. 



B. Numerical results 

We solve the set of Eqs. P^ - (|22p for a slab geometry 
and a linear distribution of the potential profile in order 
to investigate its effect on strong correlations. Although 
we do not consider long range Coulomb interactions or 
Poisson-Schrodinger coupling at this level, it is possible 
to couple the current solutions to a Poisson solver in order 
to consider more screening effects. 

The spatial distribution of the quasiparticle weights 
and the charge densities are shown in Figs. [T] and 
Fig. [2Ja)-(b) for different values of the Hubbard repul- 
sion U for a slab of width = 100. Mathematically, the 
existence of a potential profile causes charge distortion in 
the system and because of the nature of the Gutzwiller 
renormalization factors that have a minimum value at 
half filling(ni = 1.0) it is predicted that any charge frus- 
tration in the system may lead to larger quasiparticle 
weights when compared to the case without electrical 
field. 

For both C/ < C/^ and [/ > Uc (where C/c = 16t for bulk) 
the maximum quasiparticle weight is achieved in few lay- 
ers beneath the surface as is obvious from Figs. [Ija) and 
. For U < Uc akin to the zero electric field case^i the 
minimum quasiparticle weight is achieved for the surface 
sites. In contrast, for U > Uc the quasiparticle weight 
of the central atoms dramatically starts to drop to ex- 
tremely low values and creates a dead insulating region 
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FIG. 2. (a) Quasiparticle weight distribution, (b) charge dis- 
tribution and (c) effective potential for v = 2t. Notice the 
C//2 contribution is subtracted from the effective potential. 



as is indicated in Fig.[21[a). This is presented more clearly 
in Figs. [3^a)-(b) where we show the quasiparticle weight 
versus the Hubbard repulsion for three significant loca- 
tions (surface, near surface and bulk) for both v = Q and 
V = 2t. 

The formation of this dead zone together with the fact 
that we increased the value of the Hubbard repulsion 
from lower values may lead to charge being trapped near 
the surfaces of the slab because electrons are not able to 
tunnel through the bulk. This charge trapping prevents 
the system to exhibit a complete metal-insulator transi- 
tion even for values of the Hubbard repulsion larger than 
the bulk Uc- This result is contrary to the case without 
electrical field where there is a single U = Uc sA, which 
the quasiparticle weight is suppressed for the hole sys- 
tem. When there is no electrical field the quasiparticle 
weight is maximal in central parts as shown in Fig. [TJa). 
The distance over which the quasiparticle weight recov- 
ers its bulk value is on the order of 10 atoms which 
shows that the electrons which are located on the sur- 
face atoms suffering from the lack of kinetic energy (due 
to lower coordination number at the surface) are always 
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FIG. 3. Quasiparticle weight of various sites versus Hubbard 
repulsion for A'', = 100, (a.) v — and (b) v — 2t. 





0.12 




0.10 








0.08 


o 




■c 

CO 


0.06 


Q. 




"w 

CO 


0.04 






O 


0.02 




0.00 



v=2.075 
v=3.000 
v=4.000 



A 



10 20 30 40 50 60 70 80 90 100 
Position of atoms in z direction 

FIG. 4. Quasiparticle weight distribution versus v for U — 
16.0625t and Nz = 100. 



able to gain kinetic energy from the central sites with 
highest quasiparticle weight. Thus the surface quasipar- 
ticle weights will not vanish completely as long as the 
bulk quasiparticle weight is finite, although may have 
very low values.— 

To see the difference between cases with v = and 
V ^ it should be noticed that in the case in which 
electrical field is present the quasiparticle weight is max- 
imal in few layers beneath the surface. This is because 
of the fact that electrons at these locations are more in- 
tensely affected by the electrical field while they do not 
suffer from the lack of kinetic energy as do the electrons 
corresponding to cites which are exactly at the surface. 
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FIG. 5. Critical Hubbard repulsion for which the maximal 
quasiparticle weight is Z = 5.0 x 10^'' versus slab thickness 
for different electric fields. 



Therefore the distance between the maximum quasipar- 
ticle weight (as source of kinetic energy) and central sites 
is higher specially for higher sizes and as a result the cen- 
tral sites are not able to gain kinetic energy, moreover this 
sites are less affected by electrical field when one increases 
the width of the slab together with fixing potential dif- 
ference between edges to a constant value as we have 
consider in our model. Thus the metal-insulator transi- 
tion occur for sufficiently large width and some U > Uc 
for central sites before complete screening of electrical 
field as indicated in fig. Eljb) more precisely. This is very 
similar to the case of the interface of a bad metal with 
a strongly correlated material with U > Uc, where there 
is an insulating phase sufficiently far from the source of 
kinetic energy which is located near the surfacci ^^i^^ 

In Figs, mj^b) and [2jc) the spatial distribution of the 
charge densities and the effective potentials are shown 
for different values oi U > Uc- Both of these two quan- 
tities behave similarly to the quasiparticle weight. The 
charge density is maximum in the same location in which 
we have the maximum quasiparticle weight while for the 
sites with charge density near local half filling (n^ = I.O) 
we have the lowest quasiparticle weight and this is where 
the electrical field has the weakest effect. This confirms 
that the higher quasiparticle weight is due to a larger 
carrier density near the surface of the slab. The devia- 
tions of the carrier densities from half-filling correspond 
to larger electron density for sites with lower effective 
potential and hole density for sites for higher effective 
potential as shown in Fig. [2{c) . The charge frustration 
is responsible for nonzero quasiparticle weight for these 
sites near the surfaces of the system even for U > Uc- 
The regime of nonzero conductivity of the edge regions 
for U > Uc is similar to underdoped and overdoped Mott 
insulators but in this case we have an inhomogeneous 
charge distribution due to the presence of the electrical 
field. 

Fig. m shows the change of quasiparticle weight 
throughout the system when the voltage difference is in- 
creased. While the location of the maximal quasiparticle 
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FIG. 6. (a) Quasiparticle weight distribution of different sites 
for various slab thicknesses; (b) The charge density averaged 
over half of the slab for different thicknesses. Here U = 16t 
and V — 2t. 



weight slowly shifted towards the surface, its value in- 
creases with electric field. This in turn assures that the 
size of the central dead zone reduces when the voltage 
difference is increased. One should note that when mea- 
suring an I-V curve only in-plane the conductivity will 
show metallic behavior because the z-axis conductivity 
will be dominated by the bulk insulating layer. 

Fig. [5] shows the value of the Hubbard repulsion for 
which the maximum quasiparticle weight is Z = 5.0 x 
10^'^ as a function of slab thickness. This will give a 
lower bound for the critical Uf"''' in the presence of a 
perpendicular electric field. Jjf"-'' is higher for larger 
thicknesses and stronger fields v. Again this is related to 
the amount of charges localized near the surfaces. When 
U increases the quasiparticles corresponding to central 
parts drop faster in the thicker slabs as is indicated in 
Fig. [nja) and this causes more charge accumulation at 
the surfaces. This is because of the fact that the proba- 
bility of the electrons to tunnel through the central parts 
is being reduced which makes charge relaxation more dif- 
ficult. In other words by increasing the Hubbard repul- 
sion the system tries to screen the charges due to energy 
restrictions while on the other hand this increasing of U 
suppress metalic behavior of central part and thus hin- 
ders the charge relaxation. This scenario is expected to 
be even more relevant for thicker slabs. This can be bet- 
ter understood by considering the average charge accu- 
mulation in half of the slab which increases for thicker 
slabs as indicated in Fig. Hl^b). 
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U/t 

FIG. 7. The average charge density for one side of slab for 
iV;: = too. 



Fig. ([7]) shows the charge accumulation in one side of 
the slab as a function of the Hubbard repulsion for differ- 
ent strengths of the electric field. While the charge accu- 
mulation decreases for higher Hubbard repulsion due to 
screening, it instead increases for higher voltage values. 
Although we did not obtain a clear asymptotic behavior 
for maximum quasiparticle weight by increasing the slab 
width, it seems that one may gets an asymptotic solution 
for very thick slabs, in which the central sites may have 
a MIT at U very close to 16t. This may occur due to the 
very large distance of the central sites from the edges. 
There are two reasons that the central sites yield MIT 
very near 16t for very thick slabs: first of all the cen- 
tral parts have the lowest charge density deviation from 
n = 1.0 because of the shape of the potential profile, sec- 
ond because it is difficult for these sites to gain kinetic 
energy from the source of kinetic energy which is located 
in only few layers beneath the surfaces. By considering 
the latter facts an asymptotic solution may be achieved 
for extremely thick slabs. 

IV. CONCLUSIONS 

In conclusion, we described the Mott metal-insulator 
in a slab geometry in the presence of an external elec- 
trical field by calculating the site dependent quasipar- 
ticle weight. This is done by using an inhomogeneous 
Gutzwiller approximation which is exact in the limit of 
infinite dimensions. Increasing the Hubbard repulsion 
from lower values in the presence of an external electri- 
cal field leads to the formation of a dead insulating zone 
at the center of the thin film. The formation of the dead 
zone for U > 16t happens before the complete screening 
of the electrical field and therefore charge trapping oc- 
curs at the edge sites. This charge trapping causes the 
MIT to be shifted for edge sites in the presence of the 
external field. We therefore show that even though the 
central region becomes insulating at Uc, the surface lay- 
ers remain metallic but with a suppressed quasiparticle 
weight. From an experimental point of view our results 
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are relevant for transport measurements in thin films. In 
the presence of an external electric field perpendicular 
to an insulating film, one could use the surface states 
for transport since the charge transfer at the surface cre- 
ates two dimensional underdoped/overdoped regions. In 
the same time, transport perpendicular to the thin film 
is suppressed due to the dead insulating zone, thus pro- 
tecting the surface states from leakages. The electric field 
needed to create the surface states is also much lower 
than the breakdown field needed to pass current across 



the insulating zone. 



ACKNOWLEDGMENTS 

This work was supported by the Flemish Science Foun- 
dation (FWO-Vlaanderen) and the Belgian Science Pol- 
icy (lAP). L.C. acknowledges individual support from 
FWO-Vlaanderen. 



1 A. Cavalleri, Th. Dekorsy, H. H. W. Chong, J. C. KiefTer, 
and R. W. Schoenlein, Phys. Rev. B 70, 161102(R) (2004). 

^ Hyun-Tak Kim, Byung Chae, New J. Phys. 6 52 (2004). 

^ A. D. Caviglia, S. Gariglio, N. Reyren, and D. Jaccard, 
Nature 456, 624 (2008). 

* D. Ruzmetov, G. Gopalakrishnan, and J. Deng, J. Appl. 
Phys. 106, 083702 (2009). 

^ A. Husmann, J. Brooke, T. F. Rosenbaum, X. Yao, and J. 

M. Honig, Phys. Rev. Lett. 84, 2465 (2000). 
^ J. Hubbard, Proc. R. Soc. London 276, 238 (1963). 

M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963). 

* W. Metzner, D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989). 
^ W. Metzner, D. Vollhardt, Phys. Rev. B. 37, 7382 (1988). 

^° F. Gebhard, Phys. Rev. B 41, 9452 (1990). 

" E. Lieb, F. Y. Wu, Phys. Rev. Lett. 20, 1445 (1968). 

G. Kotliar, and A. E. Ruckenstein, Phys. Rev. Lett. 57, 

1362 (1986). 

A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg, 
Rev. Mod. Phys. 68, 13 (1996). 
" J. Bunemann, F. Gebhard, and R. Thul, Phys. Rev. B 67, 
075103 (2003). 

W. F. Brinkman, and T. M. Rice, Phys. Rev. B 2, 4302 
(1970). 



J. J. More, B. S. Garbow, and K. E. Hillstrom, User Guide 
for MINPACK-1, Argonne National Laboratory Report 
ANL-80-74, Argonne, 111., 1980. 

J. J. More, D. C. Sorensen, K. E. Hillstrom, and B. S. Gar- 
bow, The MINPACK Project, in Sources and Development 
of Mathematical Software, W. J. Cowell, ed., Prentice- 
Hall, pages 88-111 (1984). 

A. Ruegg, S. Pilgram, and M. Sigrist, Phys. Rev. B 75, 
195117 (2007). 

P. Fazekas, Lecture Notes on Electron Correlations and 

Magnetism, Series in Modern Condensed Matter Physics, 

Vol. 5( World Scientific, Singapore, 1999). 

M. Potthoff and W. Nohing, Phys. Rev. B 59, 2549 (1999). 

R. Nourafkan and F. Marsiglio, Phys. Rev. B 83, 155116 

(2011). 

H. J. Monkhorst and J. D. Pack, Phys. Rev. B. 12, 15 
(1976). 

G. Borghi, M. Fabrizio, and E. Tosatti, Phys. Rev. B 81, 
115134 (2010). 

R. W. Helmes, T. A. Costi, and A. Rosch, Phys. Rev. Lett. 
101, 066802 (2008). 



